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I . INTRODUCTION 


Three-point finite-difference discretization has formed the basis for 
the overwhelming majority of numerical solutions of the equations of fluid 
mechanics. For uniform meshes these procedures are typically second-order 
accurate in the mesh width h. A decrease in order of accuracy results for 
nonuniform grids. A wide variety of temporal or marching integration schemes 
have been developed and these include explicit (one step or two step methods) 
or implicit procedures. For the latter, which generally have better stabil- 
ity properties, the primary advantage of the three-point differencing is 
that the resulting algebraic matrix system is of a block-tridiagonal form; 
therefore, an efficient and well developed two-pass algorithm^^^ can be ap- 
plied to invert the matrix operator. 

Recently, a number of higher-order numerical methods have been proposed. 

The obvious extension is to five-point differencing which leads to a fourth- 
order accurate system. Unfortunately, for implicit integration the matrix 
system is pen ta-d i agona 1 and, therefore, the boundaries require special considera- 
tion. In addition, the truncation error is considerably larger than that found 

( 2 ) 

with the Spline and Hermi te methods to be discussed. Graves has proposed a 
five-point partial implicit procedure that simplifies the inversion process; 
although this method is inconsistent in the transient it can be useful for time 
asymptotic solutions. 

A second class of collocation procedures which are also fourth-order accu- 
rate for uniform meshes and which retain a 2x2 b lock-tr i d i agona 1 form for the 
governing matrix system have recently been proposed. These Hermi te or Spline 
collocation techniques treat both the functional values and certain derivatives 
as unknown at the three collocation points. These procedures generally result 
in a somewhat lower truncation error than that found with a five-point func- 
tional discretization and can be derived from appropriate Taylor series ex- 
pansions (Hermite) or polynomial interpolation (Spline). In the former category 


The blocks are 2x2 for the fourth-order methods and 3x3 for the sixth-order 
methods . 


(3) (^) 

we have the Fade approximation of Kreiss or so-called compact scheme , 
the Mehrstel lung formulation and Hermit lan finite-difference development 
of Peters^^^ In the latter group are the spline collocation methods de- 
scribed by Rubin and Graves^^^ and Rubin and Khosla^^^. In addition, a 
spl i ne-on-spl I ne techniques is shown to result from a hybrid formulation. 


The purpose of the present analysis is l) to clarify the relationship 
between the various Hermite developments, 2) to point out an apparent in- 
consistency in the Peters' formulation, 3) to derive the Hermite tridiagonal 
system for a nonuniform mesh, since all previous developments are for uniform 
meshes, k) to extend the Hermite philosophy in order to develop a variable 
mesh sixth-order tridiagonal procedure, 5) to briefly review the spline 
interpolation method, develop this collocation procedure for several new 
polynomial forms resulting in block-tridiagonal systems and to demonstrate 
that, in fact, all of the results obtained by the Hermite Taylor series de- 
velopment can be recovered by appropriate spline polynomial interpolation. 
Comparative solutions are presented for the boundary layer on a flat plate, 
boundary layers with uniform and variable mass transfer, the nonlinear 
Burgers equation, and the viscous incompress i b le Navier-Stokes equations de- 
scribing the flow in a driven cavity. Finally, 6) the polynomial Interpola- 
tion procedure is used to develop higher-order temporal integration schemes. 
Solutions are obtained for the diffusion equation describing the impulsive 
motion of a flat plate (Rayleigh problem). 
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II. POLYNOMIAL SPLINE INTERPOLATION 


Consider a mesh with nodal points Xj such that 


a = X < X- < ... < X., < X.,.- = b 
o 1 N N+1 


Define the mesh width hj = Xj - with o. = a=hj_j^^/hj. Consider a func- 
tion u(x) such that at the mesh points Xj , we specify ^ j • 

purposes of the present analysis we define the polynomial spline S(xj;n,k) = 

S(n,k), such that at the mesh points x. we prescribe S(x.;n,k) = u.. S(n,k) 

th J J 

is an n order polynomial defined on any interval [j-1,j] and in the set 

n— k 

C [a,b]; k is defined as the deficiency of the polynomial spline; i .e. , we 
are considering an n^*^ order polynomial having n-k continuous derivatives on 
[a,b] . 


The so-called simple spline^^^ has deficiency k=l . The familiar cubic 
spline is a cubic polynomial of deficiency one or S(3,l)- For a more detailed 
discussion of the properties of polynomial splines see, for example, Reference 
( 1 ). 


Cubic splines have been widely used for curve fitting and interpolation 
purposes, but only recently has spline collocation been adapted for the numeri- 
cal solution of ordinary and partial differential equat i ons ^ ^ ^ . 

These procedures have been applied to the equations of fluid mechanics by Rubin 
and Graves^^^ and Rubin and Khosla^^^. In these papers the spline collocation 
technique is described for the basic cubic spline S(3,l) as well as a higher- 
order accurate quintic spline of deficiency three S(5,3). The former has been 
termed spline 2 and the latter spline h. In addition, in Reference (8) it is 
shown that the basic three-point finite-difference discretization formulae can 
be obtained by considering the quadratic spline of zero deficiency, i.e., S(2,0), 


The general spline interpolation procedure of References (7) and (8) can be 

order polynomial is specified on the interval 


summarized as follows: an n 

The n+1 constants are related to the functional values u 


■1 


u. , 
J 


as 


well as certain spline derivatives m. m., M. M.. m., M. are the spline 

j-1 J* j-1’ J j’ J 
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derivative approximations to the functional derivatives u*(Xj), u‘*(Xj) respec- 
tively. A similar procedure Is considered on the interval Continuity 

of derivatives is then specified at Xy This process results in two coupled 

equations for mj , Mj , j=l , N. Boundary conditions are required at j=0 

and j=N+l. The system is closed by the governing differential equation for 
u(Xj), where all derivatives are replaced by t-heir spline polynomial approxi- 
mations my M j . The details of this procedure for spline 1, 2, k are given in 
References (7) and (8) where a variety of explicit, implicit, two-step, relaxa- 
tion and ADI methods are explored. 

This spline procedure can be applied to other polynomials of other orders 
and deficiencies and thereby a variety of systems can be derived. Since the 
equations of fluid mechanics are second-order we restrict our attention to 
polynomial splines defined solely by the functional values and the spline first 
and second derivatives at the nodal points. In addition, only those polynomial 
splines resulting in at most 3x3 block-tridiagonal matrix systems are considered. 
In this regard, in addition to splines 1, 2 and i.e., S(2,0), S(3,l), S(5,3), 
which have previously been described, the equations governing the spline deriva- 
tives mj , Mj for S(4,0) and S(5,l) have been evaluated; All of the governing 
systems for the various procedures are now specified. The spline polynomial 
on [j-l,j] is also discussed. Consider the polynomial spline on [j-l,j] 


n 

S (x ; n , k) = Z A. t ' (l ) 

i=o 

at the nodes specify 

S(Xj_^,n,k) == ; S(xj;n,k) = Uj (2) 


k 



Jn addition we require some or all of the following conditions: 


S'(Xj_^; n,k) = (3a) 

S' (xj ;n,k) =nij (3t>) 

S"(Xj_i;n,k) = (i,a) 

S"(x.;n,k) = M. (i,b) 

where m ^ , Mj are the spline derivative approximations of u ' (x) , u"(x) , re- 
spectively. The interior point spline equations are as follows: 

S (2 , O) or Spline 1 - Conditions (2) and either (3a) or (3b) are 
specified for the constants A. in (l). 

S(x;2,0) = Ujt + (l-t) + (uj - Uj_^ + h^mj t(l-t), (5a) 

where 

t = (x-Xj^^)/hj. 

2 2 

nij = + (a -l)uj - a Uj_^)/(o(1+a)hj) 

Mj = 2(uj^^/a - (1+a Uj +Uj_^)/(hj (l+a)) (5c) 

S(3, l) or Spline 2 and S(5,3) or Spline U - Conditions (2) and 
either (3) or {k) determine the constants A. for S(3,l), For S(5,3) we require 
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(2) and (3), but in lieu of (i|) , and in order to increase the accuracy of the 
spline second-derivative approximation, we specify 


M. = K. + G. ; M. = K. , + G. , 
J J J J-1 J--1 J-1 


(6a) 


where 


Gj = - (1+a) K. + a Kj_^)/6 


(6b) 


and 


A = (l+a^)/(a(l+a)^) 


(6c) 


In this formulation S(3,l) is recovered from S(5,3) when A = 0 so that Mj=Kj. 
The polynomial spline is given by 


s(x;5,3) = K, , (1-t)3 h?/6 + K.t^ h':/6+(u, , - K h^/6) (l-t) 

J"« J J J J * J“* J 

+ (u. - K. h?/6) t+G. t^ (l-t)^/2 + G, , t^(l-t)^/2 
J J J J J“1 


3 2, 


(7a) 


Recall that S (x . ; 5,3) S(x.; 3,1) as A 0. In addition, we have for 

both polynomial approximations 


+ 2(l+a)K^. + a = 6(uj_^^/a - (l+a ^)uj 


+ u._,)/h. 


(7b) 


m- 

J 


= h. 
J 


(Kj 


+ 0 


5Kj.,)/3 


(u. - 




(7c) 
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mj (Kj + 0.5Kj^,)/3 + (uj^, - UjJ/I'j+i 


Combining (7c) and (7d) we obtain 


a m. , + 2(l+a) m. + m = 3(u /a + u , (a~a ) 
j-1 J J+1 J+« J 


a uj.,)/hj 


As shown In References (7) the relations (7c), (7d) are generally preferable 
to (7e) as they provide a direct evaluation of and, therefore, m^ can be 
eliminated in favor of K j , u ^ . Also for nonuniform meshes where (7c), (7d) 
and (7e) for mj are third'order accurate spline representations of u ' (x) , the 
truncation error is increased if (7e) is applied. For spline 2 the governing 
matrix system is tridiagonal for M., For spline h a 2x2 bl ock- t r i d i agona 1 sys- 
tem results. 

S (^ ,0) “ The A. in (l) are evaluated from (2) and either (3) and 

(4b) or (4) and (3b). Two distinct polynomial splines are generated. We 

1 2 

designate the former as S (4,0) and the latter as S (4,0). The polynomial splines 
in each case have been obtained. The important spline equations generated from 
these polynomials are now presented. 


A : $ (4,0) 


S'(x;4,0) = Uj t + [6(uj-Uj .j)-3hj (irij+nij .j) + t 

- [8(uj-Uj_^ )-hj (5nij + + h? Mj] 


t- ^ 
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For 


M - 

IT ^ " 


2 m. +m 2 m.+l 2 , 

J J*' . / J . cr -1 \ 

+ 1 o /i4.„\u — — + m . - 0 m . , ) 

j+1 IZU+a^h a a j j-1 


3h 


(8b) 


l+ 2 a 


12 a j 


.. j 

^+1 


M. = 


2 m. +m. - 

I L± + 

3hj 12 (l+a)hj ^ a 


/ j +1 . a “1 ^ 

^ + m , - 0m . . ) 

a J j-l' 
u.-u. . 

- -J J-T 


(8c) 


The continuity of Mj leads to 

m. + (1+a)^ m. + a^m. = |_ ^ [I±l£ ^ + £zl (i+a)3 

J+1 J J-1 1+a hj a J+1 a ' ^ j 


“ a (2+a) 


a=1 it can be shown that S^(i),0) is equivalent to S(5,3). 


(8d) 


•B : S^(4,0) 


h^ 


S (x; 4,0) = Uj_^ + [2(uj-Uj_^) - hjiTij + -g-l- (Mj-Mj_^)] t 


h? 


+ t - [2(u.-u._^) - 2h. m. (M.-tM..,)] r 


2M.+M. 
L_ 

J j 6 


+ [(u.-u. .) - h, m. + — ^ h? ] t^ 


(9a) 


_ “j “j-1 I "j r1+4(j 1 ,, , 1+2a „ T 


h. 

J 


h. 

12 


(9b) 


8 


m. 

J 


= .!j±iZl 



[(4+o) 


M , + 
J 


2-fa 

1+0 


”j+i 


2 

o 

1+0 


Mj.,1 


The continuity of 


m. 

J 


leads to 


0^+0- 1 
120 


j+1 


o 2 2 

. 0^+^0 +^0+1 . 1+0-0 

TzB "• ~12 





“J * “j-i> 


(9c) 


(9d) 


It is significant that for all of the polynomial splines considered thus far 
the governing system consists of a tridiagonal form for m. or (Mj) and a pair 
of algebraic relations for the other spline derivative Mj (or ^j)- Therefore, 
at most a 2x2 b 1 ock- tri d i agona 1 system results. We shall now consider the 
simple quintic spline S(5,l). For the first time the governing interior point 
system enlarges to a 3x3 block-tri di agona 1 in terms of both of the spline de- 
rivatives mj and M ^ . The simpler algebraic relations no longer appear and, 
therefore, the final matrix system will be somewhat more complex. On the 
other hand, increased order of accuracy is achieved. 


S (3 , 1 ) - The six A. in (l) are evaluated with all of the 
conditions (2-A) . Continuity of the third and fourth derivative leads to the 
following interior point equations: 




S(x;5,l) = + hj mj^^ ^ 2 ^j -1 ^ 




+ [10(uj-Uj_^) - hj (^mj+6mj_^) + -^ t*^ 

h? 

- [15(uj-Uj_^) - hj(7nij+8mj_^) + -^ (2Mj-3Mj_^)] t 

h- 

+ [6(uj-Uj_^) - 3hj (mj+iHj ,| ) + .j)] 


(10a) 


9 


+ 8(l+a^) ntj + 7<J^ “j " 


+ ahj + I (a^-1) Mj ~ Mj_^) 


(10b) 


and 


1 +i:±2Lm. -1m. =iJ(^-J±f u. -hu. j 

20 j+1 2 o j 2 j-1 ^2 ' ^3 „3 J J-1 

j 


^ 1 r '* ^ c <J^-1 ^ J. 1 

hT ;T "’j+l ^ ~r "’j 

JO-' a 


(10c) 


Other polynomial splines can be considered, however, for polynomials of fifth 
or lower' order, the spline formulations presented herein appear to be the most 
efficient. For higher-order splines, we require that the third or higher-order 
spline derivatives be specified in the evaluation of the A. in (l). These 
formulations are not discussed here, although the tridiagonal sixth-order ac- 
curate system for Mj derivable from S(6,0) is presented later in this report. 


For the polynomial spline formulations 

T(h.) for the various spline derivatives m, 
J J 

recall that 


presented here, the truncation errors 
and Mj are depicted in Table I. We 


m. = u' (x.) + T(h.) 

J J J 

M. = u"(x.) + T(h.) 

J J J 

For completeness the truncation errors T(h.) are also given for the five point 
finite-difference discretization with a uniform grid. Note that although 
these errors are fourth-order, they are somewhat larger than those obtained 
with any of the fourth-order polynomial spline formulations. 

For 0=1 the minimum truncation errors of the fourth-order methods are ob- 
11 2 

tained with S(5,3) and S (4,0). S (4,0) and S (4,0) retain fourth-order accuracy for 
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TABLE 1 

TRUNCATION ERROR OF SPLINE DERIVATIVES 


Uniform Mesh {0=1) Nonuniform Mesh 

m. M, m. M. 

J J J J 


S( 2 , 0 ) 


k2 . 

n / 1 Vv 

12 h 

V- >j 

■■h-"- — ■ ■■" ■ " — 

' ^ 

: a-1 , f ^ ^ iV , iV4 

S( 3 , 1 ) 

/ Vx 

T 5 o 

u2 • 

h / I Vv 

11 )j 

* (a-1) , 3 / i v^ 

“ 24 ''j >j * 

2 

(1-a+a ) /,,v^ 

— ^ 

1 . 3 h! . 
1+a 12 

S( 5 , 3 ) 

/ Vv 

T 5 o 

h / vu 

Same as ( 3 , 1 ) 

7h? 

(l+a ) (0-1) (u )j 

■" ^ ^350 ■*■ To 5 o~ h 

s’(A.o) 

/ v» 

T8O 

A 

h / vu 

350 

0 ^ (l+a)^ A f V. 

^ 'j ’j 

A 2 

/ 1 \ j 3 + 5 fJ+ 3 cJ / V\ 

120 2 

1 +0+0 

^ 0^-0^+1 , vi^ 

* 2^^, 120'“ >j 

s^(A,0) 


h / VI . 

m 

3 (0^+1 ) ' 5 (1+0+0^+a^ 
720 

2 h 3 

(0-1) (2+50+20^) Y^ (u'')j 

h^» 

+ [3(a^-D^ + 0(20^-0+2)] ^ (u''')j 




Uni form Mesh (a=1 ) 


m . M . 

J J 



5 - pt. Finite-Difference: 
a=l 


Herml te 6 : a=1 ; 


*This is obtained from (7c) or 
72. 


TABLE I (Concluded) 


Non uni form Mesh 


m. 

J 


^ (H6a+6a^+0^) (u^‘) 


15(l+a^) 


3 


6 , vi i 


5^ j 




3 ^ hV 

0 ( I +CT ) J / V U 

-u 7 - 72 Q^" 


X 12h 30 ' 

I6(u,^,+u,_ -2u.) - (u,_^,-2u,+u,_ ,) ,4 


u = 

XX 


j: .... j ii^r L rjsii. + hi (^^>) 

i 2 h 2 9 ° ^ 


m. = u + h^(u'^' ' ) ./9^50 ; M. = u + h^(u^* * * ) ./ 6636 O. 

J X J J XX J 


(7d) . If evaluated from (7c) 24 In the denominator changes to 










nij even with a variable mesh. The other fourth-order polynomial splines lead 

to third-order accurate formulas for with a 1. S(5»l) is sixth-order for 

m. with a=l and fifth order wi th a / 1 . M. is fourth-order in both cases. A 
J J 

sixth order formulation for is discussed later in this paper. From Table I 

we see that even with h=1.0 there is a significant reduction in truncation 
error with the higher-order methods. This is due to smaller numerical coef- 
ficients in the error terms. 
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III.. TAYLOR SERIES FORMULATION - 
HERMITE COLLOCATION 


1, Compact Formulation - As discussed previously, higher-order 
finite-difference equations can be derived from Taylor series expansions. For 
a uniform grid, fourth-order accuracy is achieved with a five point expansion 
formula. The resulting system is pen ta-d i agona 1 with implicit integration 
procedures. Recently, Orszag and Israel i^^^ have reported an idea due to 
Kreiss in which a compact or Pads approximation transforms the penta- 
diagonal system for the functional values at the nodes to a 3^3 block-tri d 1 agonal 
system for the functional values and their derivatives at the nodes. 


It 


has been observed 


( 3 ) 


that with 


m. = 
J 



M. = (u ) 

J XX 


D 


u . 
J 




Du. = 
J 


(u.-U. , 

J J-1 



(11a) 


D°Uj = 2(Uj^j/o - (1+a) Uj/a + u^ _.|)/((l+a)hj) , (11b) 

d' U j = (Uj^^/a + (a^-1) Uj/a - aUj_^)/(l+a)hj) , (11c) 

for a uniform mesh , the five-point difference discretization is of the form 

n,j = (1 - dV) (u.^,-u._p/2h 
2 

Mj = (1 - J2 C^'*"D'uj) 

The truncation errors are given in Table I. These expressions can be re- 
written with a Fade or compact approx ima t ion such that 



(12a) 


or 


m. 

J 


1 + d'^d" 


d‘*'d‘u 


M. 

J 


1 + ^ D D 


+ - 

(l + -g- do) irij - (Uj_j_^-Uj_^)/2h 

h2 + - +- 

(1 + ^ D D ) H. = D D Uj 


(12b) 


(13a) 

(13b) 


This results in a fourth-order block-tridiagonal interior point system for 

the function u. and the deriv=^tives m., M.. As before the system is closed 
J J J 

with the differential equation and appropriate boundary conditions. 


The system (13) has appeared in a number of places and is termed com- 
pact^^^, Mehrstel lung^^^ and Hermitian^^^ differencing. The system (13) is 
fourth-order with a somewhat smaller truncation error than the five-point dif- 
ference equations. 


The E quations (13) have previously been observed in the spline analysis pre- 
sented herein. For a uniform mesh, (l3a) is equivalent to (7e) in S(3,l) and 
(l3b) corresponds to (9d) in S (^,0). Therefore, this compact formulation is 
the result of two different polynomial spline formulations for m ^ , M ^ . The 
derivation (13) does not provide the simpler expressions (7c, 7d) or (9b, 9c) re- 
lating the derivatives M.. These expressions are particularly useful in 
consideration of boundary conditions and in order to eliminate m. and, thus, 
reduce the size of the governing matrix system. Moreover, (7c) and (9d) have 
been obtained for nonuniform meshes. An extension to variable grids of the com- 
pact formulation could be quite laborious and will not be considered here. On 
the other hand, the "compact'* formulae (13) can be derived directly with a 
three-point Hermitian col.location procedure. It is the approach that will be 
discussed further. 
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2 . 


Hermitlan Collocation - Consider the finite three-point Taylor 


series expansions 


2 L 

U -.1 = u.+h,.-u + h. .u /2 + h.,-u. /6 + h..-u /2k 

J + l J J + lx J + IXX J + 1 XXX j + 1 xxxx 

u, - = u.-h.u + h^u /2 - h?u /6 + h^u /2k 
j-1 jjx jxx j XXX j xxxx 


Let 


m. = u , M. = u 
J X J XX 


(Ua) 

(Ub*) 


A. Hermi te 1 - with the operators (11) Equations (l4) can be 
rewritten in the form 

m. = (uj-Uj_^)/hj + (hj/2) [1-(hV3)D* + (hJ/12)D°]M^ (t5a) 

-"j = - (h.^/2) [1-(hj^^/3)D* - (h?^^/12)D°]Mj ( 15 b) 

These equations are identical with (9b) and (9c). Eliminating m^ the tri- 
diagonal relation (9d) is recovered. Therefore, Hermi te 1 is identical to 
2 

S (^jO). The tridiagonal Equation (9c{ represents the extension to variable 
meshes of the compact formula (I3b), since (9d) (l3b) for a = 1. Signifi- 

cantly, (9d) can be derived with a compact operator directly from three- 
point finite-difference expansion. The cumbersome Pade approximation which 
is not easily extendable to nonuniform meshes or higher order schemes is a- 
volded. This represents an alternate derivation of the polynomial spline 
equations. 


expans i ons 


B. Hermi te 2 - If we consider finite three-point Taylor series 


2 

+ h..-u + h.,.u /2 + hT^.u /6 

J + 1 J j + lxx J + 1 XXX J + 1 xxxx 


(l6a) 


m 


j-1 


m. - h . u + 
J J XX 


h?u /2 

J XXX 


h^u /6 

J XXXX 


(16b) 
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and combine with (l4) to eliminate u , (u = D^m.), we recover the Equa- 
tions (8). Therefore, Hermi te 2 is Identical with S*(^,0). The tridiagonal 
form (8d) is a nonuniform mesh extension of the “compact" formula (l3a). This 
relation is fourth-order accurate for both uniform and variable meshes, see 
Table I. 


C. Hermite 3 “ If in (16) we eliminate u and replace u by 

■ XXPv XXXX 

D^M. we obtain 
J 

+ 2(l+a)Mj + aMj_^ = (3/hj) (17) 

The block-tridiagonal system obtained with (8d) from Hermite 2 and (17) is 
termed Hermite 3- This is equivalent to what has been called spl i ne-on-spl i nes ^ ^ ^ , 
Since u is treated differently in obtaining the equations for m., M. this 

formulation does not result from a single polynomial spline interpolation. 

D. Hermite k - If we consider the b lock-tr i d i agona 1 system ob- 
tained with (8d) from Hermite 2 and (9d) from Hermite 1 for m j , M ^ , respectively, 
we have what is termed Hermite k. Once again this cannot be derived from a 
single polynomial spline and for o = 1 reduces to ( 13 ) or what has been termed 
the Fade or compact or Mehrstellung formulation. 


Hermite 6 - Uniform mesh. Consider the finite three-point 
Taylor series expansions 

2 k 

nO M a. h 'V ^ h Vi 

D Uj - Mj + ^ u + ^ u 

.2 

n^M ^ h VI 

D Mj = u + ~ u 

2 ^ 

ft* M H iv ^ h vi 

D = Mj + -g- u + ^ u 

Eliminating u*'^ and u'^' we obtain 

+ 8Mj - = 18{| (uj^^-2Uj+Uj._^)/h^ - (mj^^-mj_,)/2h} (l8a) 


From the finite Taylor series 
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D"uj - m.+-^ + -jjQ u 

V 

"’j = “xxx •" T2 “ 


D M. = u 


h2 V 


XXX 


+ -^ . u 


we obtain 

+ I6mj + 7m . = 15(Uj^,-Uj_^)/h + h(Mj^^-M^_^) 

For a variable mesh a similar procedure applies and we obtain 

2 

. :joa-l 1-1 ic 

Tg- ^-1 


1 


a .. . 38a-11-11a , . ^ 

M, , + .rr Mj 


- (■ 


36 a 

1+15a 


1 r 1+15a 

6hj a^d+ 0 ) 




I o(l5+a)x ^ o(_a+15) 1 

?(,«) '« ’ J J-'' 


1 ji|+5a 


90 ^(l+a) 


2 

7 - 

54a" 


(a-l)(l50 -1140+15) _ , a(5+4a) „ n 
~2 ""j ■" 9 0+a) "’j-l^ 


5 a j +1 5 a 5 a j 


. l^a^+lOa-10 _ 6 r(a^-a-l) ,, 

5 '"j-1~‘ h. ^ 2 j+1 

J a 


(g~l) (1+a)^(l+g^) 


h. 


+ a^(a^+a-l) " [(a+l)(g-2) - 3(a“l)(a+l) 

+ (2a-l)(a+1) Mj^^] 


The b 1 ock“ t r i d i agona 1 system (l8a), (l8b) for Mj , mj is termed Hermite 6, This 
is a sixth-order accurate formulation when combined with the differential equa- 
tion and appropriate boundary conditions. Simple relationships such as (7c) 
between m. and M. no longer exist and once again these formulae are not deriv- 
able from a single polynomial spline. (l8b) is the uniform mesh (a=l) analogue 
of the form (lOb) for m. found in S(5>1) for a variable mesh. Although we have 
not carried out the details, based on previous experience, we expect that (l8c) 
for a nonuniform mesh is derivable from some form of S(6,0). The truncation 
errors for Hermite 6 (l8a,l8b) are given in Table I. 


(18b) 


08c) 


(I8d) 
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Therefore, it is possible to derive the polynomial spline results of Sec- 
tion II with an Hermitian discretization procedure. Moreover, hybrid systems, 
which represent approximations resulting from multiple spline formulations, 
can also be conceived. One of these hybrid systems Is the variable mesh ex- 
tension of the so-called Fade or compact differencing scheme. The truncation 
errors for all possible systems can be obtained from Table I. Finally, the 
hybrid systems result In a b lock- tr i d iagonal form of m j , Mj - The simpler re- 
lations relating m, directly to M. found in the polynomial spline formulations 
are not obtained. This concept has been extended to a sixth-order system in 
Hermite 6. Higher-order approximations have not been considered. 


3 • \X _P pJ yn om i a 1 I n t e rpolation - Peters Me t hod ^ ^ ^ - Recently, 

Peters has presented an Hermitian differencing procedure for uniform meshes 

which also leads to the '’compact" or Hermite 4 3x3 bl ock-trid iagonal system 

for u j , mj and M ^ . Peters has then carried out a reduction process that ap- 
pears most attractive as a single system for Uj results; however, it can be 

shown that the resulting system is inconsistent with the differential equa- 

tion and as such results in an attendant loss of accuracy. A brief description 
fol lows • 


Consider the interpolation polynomial S(x) on the three-point interval 

[j-1, j + 1]. 

S(x) = (t^+t)/2 + Uj (1-t^) + Uj_, (t^-t)/2 


+ ot (1-t^) + (1-t^) 


as before t = (x-Xj)/h. a and 3 are two free parameters which are assumed con- 
stant on [j-1 , j + 1 ] . 

I 1 1 

With mj = S (Xj) and Mj = S (xj) we obtain 

ntj = d' uj + a/h, = hD°Uj + D”uj - 2(a+p)/h 

ri A 

m. , =-hD u. + D u. - 2(a~3)/h 

J-1 J J 


(19) 
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and 


Mj = D°Uj + 2B, = D°Uj - (6a+10B)/h^ 

= D°Uj + (6a-10e)/h^ 

A. Eliminating a and 3 from (19) we obtain 

m..i + 4tn. + m. , = 6D u. 

J+1 J J-1 J 

This is precisely (I3a) or (7e) . 

B. Eliminating a and 8 from (20) we obtain 


( 20 ) 


• (21a) 


+ 10M. + M. , = 12D°u. (21b) 

J+1 J J-1 J 

This is precisely (13b) or (9d) and, therefore, with (21a) and (21b) the com- 
pact or Hermite 4 discretization is recovered. 


Peters has not considered the tridiagonal Equations (2l) but instead has 
evaluated the differential equation under consideration at the three points 
Xj Xj , • The derivative approximations at these points are given by 

(19) and (20). This leads to three equations for the five unknowns a, 3, 

Uj , ^j+1 • ^ ^ eliminated and a single tridiagonal system for 

u. results. 

J 

However, ot and 3 can be determined independently from (19) and (20) re- 
spectively, and these results effectively imply different polynomials S(x) 
leading to the tridiagonal system (21). When a and 3 are evaluated from Peters 
tridiagonal equations the resulting values are inconsistent with those leading 
to (21a) and (21b). This can be shown in the following simple example. 


Consider the differential equation 


u + 

X 


u = 0 

XX 


( 22 ) 


Using (19) and (20) and evaluating (22) at j-1, j, j+1 foil owing Peters, we 
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I 


obtain for h=1 

(D*+D°) Uj+ 23 + a = 0 

(2D°+D*) Uj - 123 - 8a = 0 

Du.** 83 + — 0 

J 

From these equations we find 


8a = - (5D*+/»D°) Uj 

Substituting this expression into (I 9 ) we obtain for mj 

m. - (| 0 * - 1 D°) .ij 

It is clear that we do not recover even the leading term in the expansion for 
u*(x). A similar result is found for Mj . 

Numerical experiments with Burgers equation 


u^ + (u - v) u = vu 
t 2 X XX 


have shown this inconsistency- Both the viscosity (v) and convection (u - y) 
coefficients are effectively modified. For large values of v the effect of 
the inconsistency is diminished. For small v values, diagonal dominance of the 
resulting matrix system is lost. Therefore, we believe that the reliability 
of the Peters reduction procedure is questionable. 
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IV. EXAMPLES 


t. Similar Boundary Layer; Zero Mass Transfer - As a first test of the 
various polynomial spline or Hermi te formulations considered in the previous 
sections, solutions have been obtained for the similarity equations governing 
laminar boundary layer behav i or ^ ^ . 

2 

II I I ^ I 

u + fu + 3(l-f ) = 0, u = u(n), f = u 


Primes denote differentiation .with respect to n, where n = y(Re/2x)^; Re is 
the Reynolds number; y is measured normal to the surface and x along the sur- 
face. The respective velocities are v and u. We approximate the derivatives 

I II I 

u., u., f with m., M. and m., respectively, so that the governing Equations 
J J J J J 

(23) become 


M. + f. m. + = 0 

J J J J 


The additional equations for m., M., m. are given in Sections II and III for 
each of the polynomial interpolation procedures. The systems are closed with 
the boundary conditions at the surface y=0 (j=l) and the edge of the boundary 
layer y=y^ or j=N. 


The additional boundary conditions on m.j , m^, M^ , M^^ are obtained from the 
Equations (23, 24), the spline formulas typified by (8.b) or (8c), or from the 
Hermi te expansions (l4) and (15). The boundary conditions used here have 
truncation errors that parallel those for the interior systems shown on Table 
I. For spline 2 and spline 4, boundary conditions have been discussed in 
greater detail in References (7, 8); however, only third-order conditions were 
used for the spline 4 calculations so that the present results are somewhat more 
accurate. 


The results of the polynomial spline calculations are presented in Table 2 
for a variety of uniform and variable meshes. The notation a=1.5/1 means that 



TABLE 2: SIMILAR BOUNDARY LAYER SOLUTION! f"lO) 


iS=0 



hz i 

0* 

SPLINE 2 
S(3,l) 

S^4,0) 

HERMITE 4 
(COMPACT) 

SPLINE 4 
S{5,3) 

HERMITE 6 

6.0 
20.0 
1 6 .078 

1 6 .063 

0.1 
1 .0 
0.1 

0.05 

1.0 

1.0 

1.5/, 

1.8/, 

0.468634 

0.475357 

0.464325 

0.462008 

0.469597 

0.479359 

0.471666 

0.469926 

0.469600 

0.473602 

0.469600 

0.467960 

0.470025 

0.468885 

0.469600 

0.469690 


( 12 ) 

f(0) »0.469600 (ROSENHEAD ) 


)8* I 


6.0 

0.1 

1.0 

1.23227 

1.23260 

1.23258 

1 .23259 

1.23259 

20.0 

1.0 

1.0 

1.20612 

1.20863 

1.21260 

1.21863 

1.23242 

9.448 

.001 

1.8/, 

1.23604 

1.23301 

X 

I.2329S 



f“ (0) » 1.23259 (ROSENHEAd'**’) 




h. = min{h., 1), and that a=1.5 for h. < 1. The remarkable accuracy of 
J J J “ 

Hermlte 6 with the uniform mesh h=1.0 is noteworthy. It is apparent that 

significant improvements in accuracy are achieved by considering higher- 

order polynomial splines. 

2. Similar Boundary Layer: Mass Transfer - In order to carry out more 

stringent tests of the polynomial methods, boundary layers with surface mass 
transfer are considered. In this section, similarity solutions corresponding . 
to mass transfer of the type ^ ^ evaluated; in the following section, 

uniform injection and suction is specified; i.e., 'v constant so that the 
boundary layer behavior is non-similar. The subscript s denotes the surface 
values. With large injection It is possible to blow the boundary layer off 
of the surface, and with large suction the boundary layer becomes very thin 
and the shear stresses become quite large. Therefore, these boundary layer 
profiles are more difficult to approximate numerically, and provide more 
exacting tests of the spline and Hermi te collocation procedures. 

The equations governing the similar boundary layer with mass transfer are 
(23 “ 24). The only change Is in the boundary conditions (25) for f ^ , so that 
now f^ = K, where K < 0 for injection and K > 0 for suction. 

The results of these calculations are shown on Table 3 and Figures (l) 
and (2). The figures show velocity profiles for suction and injection, re- 
spectively. The flat plate Blasius profile is also included in order to em- 
phasize the extreme thinning of the boundary layer with suction and the blow- 
off obtained with injection. The polynomial solutions shown on the figures 

JL. 

are in excellent agreement with the numerically "exact” values of Emmons and 
(13) 

Leigh . These profiles are coincident with the polynomial solutions ob- 
tained with splire4 or Hermi te 6 and, therefore, are not specifically included 
on the figures. The finite-difference results are not as accurate and exhibit 
an erroneous overshoot for the suction case (Figure l). For the suction profile 
only two points lie with the boundary layer. More exact comparisons are shown 


*"Exact" here means six-decimal place accuracy. 
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TABLE 3: f"(0) -SIMILAR BOUNDARY LAYER WITH NON-UNIFORM 
MASS TRANSFER 


h 

0* 

K 

FO. 

SPLINE 4 

HERMITE 6 

EXACT 

N3/N 1 

O.l 

1.0 

0.5 

0.7394 

0.7394 


0.7394 

35/81 ! 

0. 1 

1.5 /| 

0.5 

0.7842 

0.7406 



7/21 

1.0 

1.0 

0.5 

0.7992 

0.7545 



3/21 

0.1 

1.5/,. 

10.0 

7.8903 

6.9817 


7. 1397 

2/21 

0.15 

I.5/|. 

10.0 

7.6869 

6.8703 



1/21 

0.3 

1.0 

10.0 

5.2677 

7.2178 

7.0425 


1/21 

0.1 

1.0 

-0.5 

0.2326 

0.2326 


0.2326 

48/81 

0. 1 

1.5/,. 

-0.5 

0.2317 

0.2321 



9/21 

1.0 

1.0 

-0.5 

0.2514 

0.2253 



5/21 

0.1 


- 1.2 

0.0041 

0.0046 


0.0047 

12/21 

1.0 

1.0 

- 1.2 

0.0009 

0.0045 

0.0048 


9/21 


NJ 

\J^ 










on Table 3* denotes the number of grid points within the boundary layer. 
A variety of results for uniform and nonuniform grids are presented. The 
polynomial solutions retain a high degree of accuracy for both the high 
shear suction and near separated injection cases. It is generally found 
that for equal accuracy, spline /f requires one-quarter the number of mesh 
points required in finite-difference calculations; e.g., with K= 0.5, 
f (O) = 0.739^ (81 points with finite-difference) and f (o) = 0.7392 (21 
points with spline 4). Similar behavior is found with Burgers equation and 
the cavity solutions to be discussed in the following sections. 

3- Non-Similar Boundary Layer: Uniform Mass Transfer - With uniform 

injection or suction V = constant, and the following coordinate transforma- 

.( 12 ) ^ 

1 1 on IS appl led : 

5 = Vg (xRe/2)^ ; n = y (Re/2x)^ 
til = (2x/Re)^ f(C,n); u = 
v=-tli^ = (2Rex)"^ (f + sn-nfn) 

The governing boundary layer equations become 


u + fu + C(f^ u - uu^) = 0 
nn n C n ^ 

and the boundary conditions are. 


at n=0 f=iC, u=0, and 


1 im n u 1 

The spline equations are 


iJ 


(f 




(26a) 


(26b) 


*The positive sign denotes suction and the negative sign denotes injection, 
is positive in both cases. 
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where (f_) = (f . . “ i O/A^, and wl th quasi -linear i zat ion 

^ I j * J • * » > J 

2 

(up = (2u. .u..- Uj. - u )/AC 

e ij ij u I ».j 


Iteration is used for the nonlinear term and the asterisk denotes the values 
from the previous iteration. The equation for is the same as given by (24). 
Once again the spline derivative boundary conditions are obtained from the 
governing Equation (26) and the derivative relations obtained with the poly- 
nomial interpolation of Sections II and Ml. 


wi 1 1 


( 12 ) 

For C > > 1, with suction the classical asymptotic suction profile 
be recovered, i.e. 


or 


u ^ 1 - exp(-V^ y Re) 
u 1 - exp(-2nC) 


(27) 


For injection, there has been some question^ ' as to whether the boundary 
layer will separate at a finite ^ location. This question will be addressed 
in the discussion of results which follov,fs. 


The solutions are shown on Tables 4 and 5 and Figures (3) through (6). 

With many mesh points all of the methods, including finite-difference, work 
quite well, see Figures (3a) and (5a). As the mesh size is increased the 
finite-difference solutions begin to deviate from the polynomial results. 

This is shown on Figures (3b), (3c) and (5b) but more dramatically on Figures 

II 

(4) and (6). The surface shear stress, f (^,0), obtained with the finite- 
difference method becomes very inaccurate for coarse meshes. For the suction 
calculation, the asymptotic solution (27) gives for ^ > > 1 

f"(0) = 2g 

Therefore, at ? = 1.0, f (O) ^ 2. The spline 4 results very closely approximate 
this value; these solutions are in all cases more accurate than the Hermite 4 
results. Table 4 presents the shear values for both the coarse and fine grids. 
Also shown is the velocity one grid point away from the surface. The asymptotic 
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TABLE 4: f“(0) NON-SIMILAR BOUNDARY LAYER WITH UNIFORM 
SUCTION 


c 

h 

n 

N 

F. 0. 

SPL 

NE 4 

HERMITE 4 


f"(O.C) 

u(C,A7)) 

f“(0, C) 

u(C,At}) 

f "(0, { ) 

0.09 

1.0 

1.0 

II 

0.5321 

0.6798 

0.5218 

0.5829 

0.5228 

0.5874 

0.49 




0.7776 

1 .0539 

0.7357 

1 .1769 

0.7369 

1.1860 

0.79 




0.9185 

1.3270 

0.8249 

1 .6775 

0.8235 

1 .7202 

0.95 




0.9833 

1.4580 

0.8536 

1.9751 

0.8491 

2.0599 

1.0 




1 .0022 

1.4970 

0.8600 

2.0745 

0.8544 

2.1777 

0.09 

0.1 

1*5 

21 

0.0628 

0.6340 

0.0576 

0.5817 

0.0577 

0.5823 

0.49 




0. 1 253 

I.3II9 

0.1 119 

1 . 1748 

0.1120 

1.1762 

0.79 




0.1761 

1.8898 

0.1556 

1.6822 

0.1557 

1 

1 .6828 

0.95 

i 



0.2038 

2.2151 

0.1792 

1.9678 

0.1792 

1.9675 

1.0 




0.2125 

2.3184 

0.1866 

2.0587 

0.1866 

;2.058l 

0.09 

0.1 

1.0 

61 

0.0575| 

0.5807 

0.0575 

0.5807 

0.0575 

10.5807 

0.49 




0.0979 

1.0167 

0.1122 

I.I78I 

0.1122 

LI780 

0.79 




0.1566! 

1.6804 

0. 1563 

1.6902 

0.1563 

1.6900 

0.95 

1 



0.1806 

1.9629 

0. 1802 

1.9790 

0. 1817 

1.9970 

1.0 




0.1882 

2.0526 

0. 1877 

j 2. 0709 

jO. 1877 

j 2.0707 
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TABLE 5: f"(0) NON-SIMILAR BOUNDARY LAYER WITH UNIFORM 
BLOWING 


c 

h 

<r 

N 

F. 0. 

I HERMITE 4 

SPLINE 4 ' 


i ■! 



u(f,Ai7)|f“(0,$) 

u(f iM) 


u(f, A17) 






N> 


N =61, h=O.I, (J = 1.0 


I 

I 

F.D., SPLINE 4, HERMITE 4 




















solution (27) gives at S = 1.0 


u(0.l) = 0.1812 
or 

u ( 1 . O) = 0.8647 

Once again the spline 4 results are best. 

Detailed injection solutions are shown on Table 5* For the very fine grid, 
there was no indication of separation as inferred in Reference (l4). This 
was true for all claculat ions . The shear decreased but never vanished. For the 
coarser grid the finite-difference solution did lead to a separation point, 
but the polynomial solutions still did not separate. This behavior is also 
depicted on Figure (6b), The conclusion obtained from these results would 
appear to be that separation does not occur with uniform injection; instead, 
the shear decreases asymptotically to zero for large values of 

4. Burgers Equation - Since many results have been obtained in earlier 
studies of the nonlinear Burgers equation, only a brief discussion of 
the higher-order solutions is presented here. The governing differential and 
spline equations are, respectively, 

+ (u - 0.5) u = vu 

t ' X XX 

and 

(u^^^ - u*?)/At + (u? - 0.5) m. = vM. (28) 

J J J J J 

The boundary conditions are 


= 1.0 


u,=0.0 


The boundary conditions on m ^ are obtained from (28) and/or the spline 
derivative relations. More details on these boundary conditions are given in 
References (4), (7) and (8). 


Typical solutions are shown on Tables 6 and 7 for v = 1/8 and v = 1/16, 
respectively. In both cases the fourth-order methods represent an improvement 
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TABLE 6: 

SOLUTION OF BURGERS EQUATION, vl/8, <r« 
31 EQUALLY SPACED POINTS 



^ FO. 

SPLINE 2 

HERMITE 4 SPLINE 4 

s' 

EXACT ' 

0 

0.5000 

0.5000 

0.5000 

0.5000 

0.5000 

1 

CM 

• 

0 

1 

0.6999 

0.6860 

0.6903 

0.6900 

0.6900 

• 

0 

1 

:0.8447 

0.8290 

0.8322 

0.8322 

0.8320 

1 -0.6 

i 0.9269 

0.9160 

0.9169 

0.9170 

0.9170 

! -0.8 

0.9673 

0.9620 

0.9608 

0.9609 

0.9610 

1 

-1.0 

0.9857 

0.9830 

1 

0.9820 

1 

0.9820 

0.9820 

-1.2 

0.9938 

0.9930 

1 0.9918 

0.9918 

0.9920 

- 1.4 

0.9973 

0.9970 

0.9963 

0.9963 

0.9960 

- 1.6 

0.9988 

0.9990 

0.9983 

0.9983 

0.9980 

- 1.8 

0.9995 

0.9990 

0.9993 

0.9993 

0.9990 




TABLE T. 

SOLUTION OF BURGERS EQUATION, v • 1/16, <r • 1.0 
19 EQUALL/ SPACED POINTS 


\U 

Fa 

SPLINE 2 

HERMITE 4 

SPUNE 4 

EXACT 

0 

0.5000 

0.5000 

0.5000 

0.5000 

0.5000 

0.2 

0.9000 

0.8231 

0.8383 

0.8356 

0.8320 

0.4 

0.9878 

0.9641 

0.9593 

0.9617 

0.9608 

0.6 

0.9986 

0.9952 

0.9922 

0.9916 

0.9918 

0.8 

0.9998 

0.9995 

0.9981 

0.9982 

0.9983 

1.0 

1.0 

0.9999 

0.9997 

0.9996 

0.9997 

1.2 

1.0 

1.0 

0.9999 

0.9999 

0.9999 









over the second-order finite-difference and mixed-order spline 2 procedures; 
however, as with the previous examples, the spline k solutions are most ac- 
curate. Similar results were obtained with nonuniform grids^^^. As a general 
rule of thumb, it was found that the spline 4 solutions required one-quarter the 
number of mesh points of the finite-difference calculations in order to achieve 
equal accuracy; (e.g., with v = 1/8 and x = - 0.4, u = 0.8351 with 101 points 
(finite-difference) and u = 0.835& with 27 poincs (spline 4))l. This fact was 
discussed earlier for the similar suction boundary layer and will be demonstrated 
in more detail in the following section describing the flow in a driven cavity. 

5. Incompressible Flow in a Cavity - As a final test problem the laminar 
incompressible flow in a driven cavity is considered. This problem has been 
studied extensively by many investigators, see Reference (15). The governing 
equations in terms of a vorticity and stream-function system are: 

(29) 

(30) 

where \p is the stream function, ^ is the vorticity; u = and v = - i|; are the 

y X 

velocities in the x- and y-d i rect ions , respectively. The boundary conditions 
and geometry are shown in Figure (7). 


XX 


+ 


yy 


= c 


+ uC + VC = *5— 
t X y Re 


+ C ) 

XX yy' 


Solutions of (30) are obtained with a pred i c tor-corrector procedure describ- 
ed in References (16)- (17). For the Poisson Equation (29) a modified version 

( 1 8 ) 

of Buneman*s direct solver' is used. The spline approximations of (29, 30) 
in non-divergence form are 


fii' 


"tj- 


C. 


At 


^i j 


(m^ 


n+1 


U 


+ V. 


I J 




n+1 


Re 




n+1 


( 31 ) 


where J!.|j and L.j denote the polynomial approximations of ( )^ and ( 


1*5 







respectively. The superscripts denote the specific function. The spline 
boundary conditions are obtained by satisfying (31) at the walls. The de- 
tails of this procedure have been previously described in References (^) , (7) 
and (8), 

Solutions have been obtained for a square, cavi ty with Re = 100. A 17x17 
point mesh has been considered. The results are shown in Tables 8 through 12 
and Figure (8). Table 8 compares the Hermite A and spline 4 solutions with 
results obtained in previous studies^^^. The maximum value of the stream func- 
tion, the corner point velocity u and the vorticity at the mid-point of the 
upper moving wall are presented. It is significant that with a 17x17 mesh 
the spline 4 solutions parallel those obtained with a 57x57 point finite- 
difference d i scret izat ion. Moreover, the spline 4 results are very close to 
the extrapolated values as projected for the non-divergence equations. This 
behavior carries over to the velocity profile shown on Figure (8). The com- 
plete computer solutions are given on Tables (9) through (12). 





- 0.2 


NON DIVERGENCE FORM 

ED., 15 X 15 POINTS 

SPLINE 4, 17 X 17 POINTS 

F. D.. 57 X 57 POINTS 


FIG. 8 COMPARISON OF CALCULATED VELOCITY u THROUGH 
POINT OF MAXIMUM f FOR Re*IOO 





TABLE e: COMPARISON OF RESULTS FOR THE SQUARE 
CAVITY, Re =100 


C.'iiCcufuf iiiii 
Method 

tAimfHir 
of Points 

Mov itiK 

I'uint 

Vorticity 

Spline 2 

iS X 15 

7.1376 

Spline 2 

29 X 29 

6.6876 

Extrapolated Spline 


6. 5376 

Finite Difference 

15 X IS 

8.9160 

Finite Difference 

57 X 57 

6.6960 

Extrapolated 



Finite Difference 


6.5480 

Finite Difference • 

17 X 17 

7.3755 

Finite Difference • 

65 X 65 

6.6091 

HERMITE 4 

I5x 15 

6.927 

SPLINE 4 

I7x 17 

6.6104 

I * Divergence Form 


(a) Vorticity at the center of the moving surface. 


Calculation 

Number 

Maximum 

Method 

of Points 

Stream Function 

Spline 2 

15 X IS 

-.10529 

Spline 2 

29 X 29 

-.10432 

Extrapolated Spline 


-.10599 

Finite Difference 

IS X IS 

-.08742 

Finite Difference 

57 X 57 

-.10128 

Extrapolated 



Finite Difference 


-.10220 

Finite Difference • 

17 X 17 

-.09867 

Finite Difference • 

65 X 65 

-.10318 

HERMITE 4 

15 X 15 

- 1014 

SPLINE 4 

17 X 17 

-.102 3 

1 • Divergence Form 


(b) Maximum stream function. 


Calculation 

Method 

Number 
of Points 

u Velocity 
at Comer* 

Spline 2 

IS X IS 

-.13230 

Spline 2 

29 X 29 

-.10036 

Extrapolated Spline 
Finite Difference 

IS X IS 

-.08971 

.05730 

Finite Difference 

57 X 57 

-.06515 

Extrapolated 
Finite Difference 


-.074.38 

Finite Difference 



(Interpolated) 

17 X 17 

.02079 

Finite J)ifferenco ** 



(frtteriHilaCeii] 

6.5 X 65 

- . 07560 

HERMITE 4 

15 X 15 

-.118 

SPLINE 4 

17 X 17 

-.086 

*• Divergence Fern 




(c) Corner point u velocity 
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TABLE 9 

NON- DIVERGENCE FORM SPLINE 4 CALCULATED STREAM FUNCT40N 
FOR 17 X 17 POINTS . Re * 100 



































TABLE 10 


NON -DIVERGENCE FORM SPLINE 4 CALCULATED U-VELOCITY 
FOR 17 X 17 POINTS, Re »I00 


Y 

.0625 

.1250 

.1815 

.2500 

.3125 

.3750 

.4375 

.5000 

.5625 

,6250 

,6875 

.7500 

.8125 

.8750 

.9375 

0 

.onoo 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

.0000 

" .0000 

.0000 

■ .0000 

.0000 

.0000 

.0000 

.0625 

-.0010 

-.00 S 2 

-.0119 

-.0198 

-.0277 

-.0344 

: -.0389 

-.0905 

-.0384 

-.0330 

-.0250 

-.0159 

-.0078 

-.0022 

.0000 

.1250 



-.02^7 

-.0387 

-.052 h 

-.0639 

! -.0717 

-.0743 

-.0710 

-.0621 

-.0987 

-.0332 

-.0186 

-.0075 

-.0014 

. 1815 

-.0055 

-.01 P 9 

-.0367 

-.0'^ b 2 

-.0749 

-.0907 

-.lUlb 

-.1060 

-.1026 

-.0916 

-.0790 

-.0426 

-.0312 

-.0139 

-.0032 

.2500 


-.02^9 

-.0974 

-.0716 

-.096 U 

-.1151 

-.1300 

-.1373 

-. 1350 

-.1240 

-.1031 

-.0758 

-.0467 

-.0226 

-.0053 

.3125 

-.OOPP 


-.0559 


-.1118 

-.1364 

.-.1549 

-.1678 

-.1698 

-.1593 

-.1373 

-.104*4 

-.0665 

-.0318 

-.0080 

.3750 

-.oo^^v 

-.032 A 

-. Obl 4 

-.0027 

-.1236 

-.152? 

-.1765 

-.19^0 

-.2017 

-.1963 

-.1752 

-.1387 

-.0919 

-.0454 

-.0118 

.4375 

-.0100 


-.0634 

-.0967 

-.12 H 1 

-.1591 

-.1670 

-.2098 

-.2244 

-.2265 

-.2112 

-.1756 

-.1223 

-.0633 

-.0170 

.6000 

-.0096 

-.0323 

-.0613 

-.09?4, 

-.1235 I 

-.1537 j 

-.1823 

-.2081 

-.2287 

-.2399 

-.2352 , 

-.20/6 

-.1593 

-.0848 

-.0238 

.5625 

-.OOPH 


! -.0 S 68 

-.0^28 ! 

-.1090 1 

-.13 m 3 j 

-.1590 

-.1834 

-.2065 

-.2256 

-.2345 

-.2229 

-.1803 

-.1077 

-.0321 

.6250 

-. OOP ? 

-.027? 

-.04 H 6 

-.0484 

-.084/ ' 

-.1019 1 

-.1171 

-.1344 

-.1446 

-.1771 

-.1981 j 

-.2079 

-.1889 

-.1274 

-.0415 

.6875 

-.00 P 9 

-.027? 

-.0^27 

-.0 S 21 : 

-.0562 ' 

-.0577 

-.0594 

- . 

-.0757 

-.0950 

-.1223 ! 

-.1520 

-. 1666 

-.1363 

-.0521 

,7500 

-.0155 

-.03^3 

-.0416 

-.03 b 2 

-.0224 

-.0050 

.0114 

.02.30 

.0259 

.0163 

-.0093 

-.0512 

-.1002 

-.1219 

-.06 r »9 

.8125 

-.027 h 

-.osio 

- • 04*48 

-.OlbS 

.0233 j 

• U 6>* 1 

.1032 

.13*47 

.1449 

.1491 

.1415 1 

.0964 

.0196 

-.0620 

-.08**2 

.8750 

-.0609 

-.0^90 

; -.022 b 

.046? 

.1227 ! 

.19?7 

.254 1 

. 3039 

.3 J 89 

.3550 

.3456 ! 

.3 u 36 

.2125 

.0845 

-.06 m 8 

,9375 

-.1150 

.0(^11 

1 .2149 

.3123 

-42^3 j 

.9972 

.5546 

.4979 

.6261 

,6364 

.6248 i 

.5824 

.9992 

.3291 

.0563 

l.O 

1.0000 

1 .0000 

1.0000 

1 • 0 rj G 0 

1.0000 ; 

1 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 ! 

1 

1.0000 

1.0000 

1 

1.0000 

l.OOOO 




Kn 

ro 


TABLE II 


NON -DIVERGENCE FORM SPLINE 4 CALCULATED V-VELOCITY 
FOR 17 X 17 POINTS, Re «IOO 




























TABLE 12 

NON -DIVERGENCE FORM SPLINE 4 CALCULATED VORTICITY 
FOR 17 X 17 POINTS, Re »IOO 




.1250 

.1815 

.2500 

.3125 

.3750 

.437 5 

.5000 

.5625 

.6250 

.6875 

,7500 

.8125 

.8750 

.9375 

0 

.005b 


-.17.37 

-. 319H 

-.6755 

-.ha38 

-.7017 

-.7310 

-.6990 

-.5850 

-.6250 

-.2382 

-.0639 

.0133 

1 .0216 

.0625 


-.1315 

-.2JC6 

-.3135 

-.6061 

-.6899 

- • 5366 

-.55*.! 

-.5217 

-.6658 

-.3723 

-.2783 

-.1779 

-.1028 

-.0313 

.1250 


-.1^39 

-.3J7S 

-.3*^51 

- • 36ti J 

-.3777 

-.6088 

-.6070 

-.6057 

-.37P6 

-.3566 

-.3109 

-.2669 

-.1906 

-.1046 

1 . 181 5 

-.2375 

?635 

-.3351 

-.3551 

-.3502 

-.2783 

-.2740 

-.2V()rt 

-.3959 

-.3285 

-.3471 

-.3722 

-.3569 

-.3090 

-.1976 

.2500 




-.loan 

-.1759 

-.1668 

-.1638 

-.l6l3 

-.1898 

-.2691 

-.3497 

-.6275 

-.6852 

-.6576 

-.3583 

.3125 

-,SA?C 

-.36^7 

-.3216 

-.1359 

-.0Sfi9 

-.0038 

, 06o8 

.0371 

-.0016 

-.1256 

-,2«13 

-.6775 

-.5166 

-.6736 

-.5792 

.3750 

-,7395 

-.6??3 

-.2112 

-.0389 

,07^5 

.?0P0 

.2714 

.3229 

.2763 

,1695 

-.1082 

-.4205 

-.7636 

-.9226 

-.9088 

.4375 


-.5012 

-.1740 

.0672 

.2537 

.61-5 

,5560 

.579S 

.7152 

.5853 

.2956 

-.1995 

-.7533 

-1.20f>7 

-1.3367 


-1.1460 

-.S6?3 

-.IPOO 


.^165 

,5806 

.9088 

I.13?b 

1.2669 

1.2297 

.9319 

.3298 

-.5639 

-1.42.17 

-1.8989 


-1 .3333 


-.1270 

.366^ 

.5001 

.90H0 

1.2533 

1.5558 

1 .8320 

1.9399 

1.7966 

1.1685 

-.0003 

-1 .4*i2fc 

-2.5515 

.6250 

-1.536^ 

-.71^3 

-. 1 2nh 

.J623 

.7??H 

l.l?3b 

1 .5051 

1.9^S2 

2.2931 

2.6021 

2.6666 

2.2633 

.9822 

-1,1682 

-3.2686 

.6875 

-1.72?6 


-.1045 

.39b0 

.8515 

1.2579 1 

1.7053 

2.1605 1 

2.6123 

3.0316 

3.3650 

3.3212 

2.3520 

-.3030 

-3.8875 

.75 00 

-1.9(.6P. 


" . f(S34 

.5Pb9 

.9556 

1 .6387 

1 .8665 

2.3362 

2.8008 

3.3066 

3.7703 

6.1207 

3.7937 

1.3153 

-6,1298 

.8125 

-2.1339 

-.PO^O 

.3227 

.9688 

2 .6579 

1.900U 

2.3353 

2.7609 

3.1831 

3.6175 

6.1003 

6.5857 

6.9630 

3.6245 

-3.1817 

.8750 

-1 .<>'•00 

.7**0? 

1.9/33 

2.7138 

3.1592 

3 ,6997 

3.72M3 

3.9656 1 

6.1378 j 

6 * 3667 

6. 6368 

5.0262 

5.7603 

6.1700 

.5152 

.9375 

3.P2P9 

7.7359 

8.2772 

8.0023 

7 • 6 55b 

6.M89b 1 

6.38h7 

S.9S22 

S.6686 

5.6696 

5.6888 

5.7386 

6.7061 

8.2404 

8.2117 

i.O 

'•0.02S? 

37.1963 

18.71m/ 

),3.b5n 


b.5616 ' 

7.6909 

5.5)06 

5.1558 

5.2505 

6.8965 

8.5219 

10.9989 

17.9490 

28.9730 





























V. HIGHER-ORDER TEMPORAL INTEGRATION 


From the preceding sections It is evident that higher-order spatial ac- 
curacy can be achieved by using polynomial interpolation. However, in these 
calculations, the temporal accuracy is quite low; either first-order for a 
fully implicit scheme or second-order for the .Crank-NI chol son method. In 
order to achieve over all fourth or higher-order accuracy, the temporal ac- 
curacy should also be improved. Various multi-step methods in conjunction 
with 5 point finite-differences for space dimensions have been proposed in 

/ 1 Q on) 

the literature . These schemes are explicit and have severe stability 

restriction on the allowable time-step. 

Some of these explicit methods can be directly extended for the present 

spline collocation procedures. However, there are many drawbacks of such 

higher-order time integration schemes when applied in conjunction with the 

spatial spline procedures; e.g., (l) the stability restriction on the allowable 

time-step is even more restrictive than with finite-differences; (2) due to 

the large number of spline curve fits required at each time step, the overall 

computational times increase. For example, for a fourth-order Runge-Kutta 

method, at least four curve fits of each function per time-step must be 

evaluated and stored. In view of these restrictions only implicit higher- 

order methods are investigated in the present study. Some of the resulting 

methods can be found in the literature on numerical studies of ordinary dif- 

(2t) 

ferential equations . One such formulation has been recently described by 

(2 2 ) 

Watanabe and Flood and leads to an overall fourth-order unconditionally 
stable scheme. For the following analysis polynomial interpolation is used to 
derive several time integration procedures of third and fourth-order accuracy. 
Although no attempt has been made here to extend the procedure to higher-orders, 
the method of derivation is quite general. The utility of these higher-order 
temporal methods is counterbalanced by their sensitivity to initial conditions, 
possible ’’stability" limitations and increased computer time and storage. 

A. Polynomial Interpolation 


0 Crank-Nicholson Scheme - The following interpolation polynomial S(t), 
on the interval t._^ ^ ^ ^ satisfies the three conditions given earlier as 


(2) and (3b). 
5h 


S(t) = U.T + U._^ (1-t) + (u.-Uj_^-At fUj) T (1 -t) 


where 


/3ux 

m, = (-^) 


At = t. - t. - 
I 1-1 


Differentiating, we obtain 


^ (t=0) = 2(u. - 


u. - u, - . 

' At ' "' ° 2 

This is the weil-known Crank-N i chol son scheme which is second-order accurate 

(23) 

and is unconditionally stable 

12) Third-Order Accurate Method - The following cubic polynomial satisfies 
(2) and (3): 

S(t) = u._^ + T ni. 1 + [At (m. + “ 2(u. - u._^)] 

+ [3(Uj - (m.+2m._^)] 

As specified, the polynomial is completely determined by the values of the func- 
tion and its first derivative at two points i-1 and i; considering (32) on 

fh-r h.ii , we obtain at t.^^ 

~ 5u. - + 2At m. - - ku, + 4At m. 

1+1 1-1 1-1 I j 

m.^- = 12u. - + 5At m. - - 12u. + At m. 

1+1 1-1 1-1 I I 

Eliminating we get 


^±1^ = iV (5m,,, .8m, -m..,) 



From (33) other procedures can be developed; e.g. , if we eliminate u. , we 
obtain 


u. 


i+1 


u. 


1-1 


At 


= i (m..^ + ifm. + 




( 35 ) 


Equations (34) and (35) are typical of the Adams-Bashforth methods for time 
integration (see Reference 2T ) . It should be noted that both of these pro- 
cedures require one more storage location than the Crank-Nichol son method. 

The stability of these and other such schemes is discussed in Reference (2l). 
Unlike the C-N method, (34) and (35) are conditionally stable. The temporal 
accuracy of both is O(At^). The coefficients In the polynomial (32) can al- 
§o be evaluated to include the values at three node points (i“l, i, i+l). 

The results are unchanged. 


(36) 

(37) 


3) Fourth-Order Method - If the cubic polynomial (32) is considered 
on [i-1, i] and we evaluate the function and the first derivative at the mid- 
point t._^ = ^^i-1 find that 


u. - u. ^ 

n == 7T (m. 1 + 4m. ^ + m.) 

At 6 1-1 i-i 


and 


u. + u. 


u. 


i -1 At 






Equations (36) and (37) constitute a two-step procedure, which has recently 

( 22 ) 

been suggested by Watanabe and Flood . Once again the increase in storage 
is minimal and only one additional curve fit is required at the half time step. 
The method has been shown to be unconditionally stable and of fourth-order 
accuracy (see Reference 22). Unfortunately, there are convergence restrictions 
associated with the required iterative method of solution. 


For the simplest iterative procedure, where m._jL is treated explicitly, it 

^ ^ 2 

has been shown in Reference (22) that the solution converges only for At = 0(Ax ). 
This is quite restrictive. By utilizing specialized iterative procedures, the 
method has been shown to converge for At = 0(Ax), In the present 1 nvesti gation 
a novel iterative procedure has been developed from the polynomial interpolation 


56 


formula (32). If we add and subtract 2(m. + to (36) and treat 

tm._, ,^(m. + m._-)/ 2 ] explicitly, the procedure also converges for At = 0 (Ax). 
The use of this factor is suggested by the following relation, obtained from 
(32): 

m. + m. - . 

"i-i = " T “ ”i-i^ 

where 


The present iterative procedure is equivalent to an explicit treatment of the 
second derivatives in ( 38 ) and, therefore, has improved convergence and 
'•stability'* properties. 

B . fluadrature Methods 

An alternative approach to the derivation of these methods is through the 
use of quadratures. Let 

= f(x, u, Vu, V u) 


represent a general partial differential 
[ t. , t . ] , we obtain 

1 


u. - u. . = At 

I 1-1 



f dr 


equat i on . 


Integrating over the interval 


Polynomial interpolation is now used for f. 
1) Linear: 


f=f. T+f._^ (1 -t) 

Equation (39) becomes 


u. 

I 


“i-1 


At 


f. + f. 1 

I 1-1 


which is the C-N formula. 



2) Quadratic - f = f. t + (1-t) + [f. - - At (f^) ]t 0-t) 

where (f ) is evaluated from 3 point extrapolation as: 
t • 


\ 

(39) now becomes 


- '‘^ 1-1 ^ i -2 

2At 


u. - u. - - 

I 1-1 1 


At 


12 


(Sf, + 8f._, - f,_2> 


which is the Adams-Bashforth formula. 


J ) Quadratic: 


f = f. T + f._^ (1-t) + C t(1-t) 


('* 0 ) 


Therefore, 


At mf. = (f^) = f. - f - C 

Atmfj_^ = (f^) - ^i-1 where 

i -1 

C is given by f + f 

2 ^ ? 

W'ith (^0), Equation (39) becomes 


u. - u, - - 


and 


f. 


f. + f . , 

I 1-1 


i-i 


^ (mf. - mf._^) 


(i(1a) 

(i»1b) 


Since f. = (u.) = m. and mf. - (u ) ^ = M. , Equation (41-a) is exactly the 

I t j I I tt ] I 

one given in (38). Also (Al-‘b) can be written as: 


8 -r 


u. + u. 


or 


u. 


u. + u. - 
I 1-1 At 


i-i 


2— "T('i -^-1^’ 

- Vi> 


(42) 
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(^1-a) and (^2) are the fourth-order expressions given by Watanabe and Flood 
and derived earlier. 

As model problems we have considered the diffusion equation describing the 
impulsive motion of a flat plate (Rayleigh) and the nonlinear Burgers equation 
discussed earlier. Transient and large time solutions have been obtained. The 
Rayleigh problem is described by the following equations 

Uj. = vUyy , u = u(t,y) 

with boundary conditions 

u(t,0) = 1 , u(t,“) = 0 

and initial condition 

u(0,y) = 0 

The exact solution is 


u(t,y) = erfc — - — 

2/vT 

Typical solutions for the Rayleigh problem are presented in Tables 13 through 
17 . Due to the large gradients found for small times, the higher-order methods 
do not lead to a significant improvement in accuracy; in fact, at certain times, 
the lower-order solutions are in better agreement with the analytic values. If 
the exact values are used as initial conditions at t = 1.0, the large gradients 
are avoided and the results are shown on Tables 1A, 15 and I 6 . For small At, it is 
difficult to discern any difference with the various methods. For larger At 
values, some improvement can be observed; however, as At is further increased 
the higher-order methods exhibit some form of instability. For the third-order 
method this is apparently the numerical instability discussed earlier. For 
the fourth-order method, the difficulty Is associated with the iterative method 
of solution. Finally, on Table 17 large time solutions are shown. No signifi- 
cant conclusion can be drawn with respect to the advantage of the higher-order 
methods. For the Burgers equation similar behavior was found; moreover, there 
was no significant reduction in the number of steps required to attain the 
steady state. Generally speaking, for the examples considered here, there was 

59 



ON 
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TABLE 13: RAYLEIGH PROBLEM t 


t 

y 

_ . . . ... 

C -N 


0.2 

0. 

,10000E+0l 


.lOOOOE+OO- 

,«9a80E+00 



.20000E+00 

,lA65bE+00 



.30000E+00- 

.33927E-01 



,uoooot+oo 

,76978E-02 



,b0000t+00 

, lb7b9E-02 



,60000t+00 

,32bAiE-03 



,70000t+00- 

,63ab9E-0^1 



♦«oooot+oo 

, 1248lE-oa 



.90000E+00 

,23/llE-05 



,ioooot+oi 

,ab282E-06 


0.3 

0, 

,10000E+01 


.lOOOOEfOO 

.58130E+00 



.20000E+00 

,27l56E+00 



..30000E+00 

,92fa67E-01 



.aooooe+oo 

,28355E-01 



,bOOOOE+00 

,686b9L-02 



.60000E+00 

. l6a73E-02 



.70000E+00 

.380288-03 



.eooooEtoo 

.83881E-04 



.90000E+00. 

.181008-04 



♦lOOOOE+01 

,379728-Ob 





0, At =0.1 


3rc3 ORDER 

EXACT 

.lOOOOE+Ol 

,10000E+0l 

.53552E+00 

,527098+00 

.154768+00 

,205908+00 

.3323AE-01 

.57780E-01 

.72006E-02 

,114128-01 

,136098-02 

,156548-02- 

.26603E-03 

,143068-03 

,480058-04 

,967078-05 

.90825E-05 

t aa236E-06 

,160628-05 

.14596E-07 

.29281E-06 

,363368-09 

.lOOOOE+Ol 

,100008+01 

.57145E+00 

,605588 + 00. 

,294428+00 

,301708+00 

,994808-01 

,12134E + 00- 

,264158-01 

,5a667E-01 

.64722E-02 

,982318-02 

,141988-02 

,194578-02 

.304638-03 

,300838-03 

.'61232E-04 

,362958-04 

,122368-04 

.343308-05 

.23418f-05 

,257458-06 









TABLE 14; RAYLEIGH PROBLEM to'LO, At = Q.4 


i 

.» i 

y 

1st ORDER ' 

4th ORDER 

EXACT 1 

jO.l 

.0 

i.oonooE^oo 

l.OOOOOE+OO 

l.OOOOOE^OO ; 

l.OOOOOE-01 

7.86842E-01 

7.87405E-01 

7.87406E-01 ; 


2.00000E-01 

5.88706E-01 

5.89641E-01 

S.89638E-01 

1 ■ 

3.00000E-01 

4.17469E-01 

4.18494E-01 

4.18492E-01 

t 

1 

4.00000E-01 

2.7986AE-01 

2,80714E-01 

2.80713E-01 

1 

1 

5.00000E-01 

1.77012E-01 

1.77530E-01 

1.77530E-01 


6.00000E-01 

1.05478E-01 

1.05645E-01 

1.05645E-01 


7.00000E-01 

5.91578E-02 

5.90573E-02 

5.9058SE-02 

1 

8.00000E-01 

3.12138E-02 

3.09706E-02 

3.09716E-02 


9.00000E-01 

1.54937E-02 

1.52186E-02 

1.52192E-02 


1.00000E*00 

7,23fl3AE-03 

7.00062E-03 

7.00069E-03 

0.4 

.0 

l.OOOOOE+00 

l.OOOOOE^OO 

l.OOOOOE+00 

l.OOOOOE-01 

8.09765E-01 

8.11071E-01 

8.U070E-01 


2.00000E-01 

6.30337E-01 

6.32S89E-01 

6,32585E>01 


3.00000E-01 

4.70661E-01 

4.73294E-01 

4.73290E-01 


4.00000E-01 

3.36556E-01 

3.38983E-01 

3.38980E-01 


5.00000E-01 j 

2,30208E-01 j 

2.31999E-01 1 

2.31998E-01 


6.00000E-01 ' 

1.50519E-01 

1.51493E-01 

1.51494E-01 


7.00000E-01 

9.40487E-C2 

9.42623E-02 

9'.42643E-02 


8.00000E-01 

5.61632E-02 

i 5.58269E-02 

5.58295E-02 


9.00000E-01 

3.20700E-02 

3.14417E-02 

3.14439E-02 


l.OOOOOE+00 

1.752A7E-02 

1.68261E-02 

1.68274E-02 

0.6 

.0 

I.OOOOOE-^00 

l.OOOOOE^OO 

l.OOOOOE^OO 


l.OOOOOE-01 

8.21629E-01 

8.23064E-01 

0.23O63E-O1 


2.00000E-01 

6.52203E-01 

6.54724E-01 

6.54721E-01 


3.00000E-01 

4.99300E-01 

5.02339E-01 

5.02335E-01 


4.00000E-01 

3,681b6E-0l 

3.71097E-01 

3.71094E-01 


5.00000E-01 

2.61207E-01 

2.63555E-01 

2.63552E-01 


6.00000E-01 

1.78226E-01 

1.79713E-01 

1.79712E-01 


7.00000E-01 

1.16920E-01 

1.17523E-01 

1.17525E-01 


8.00000E-01 

7.37519E-02 

7.36359E-02 

7.36384E-02 


9.00000E-01 

4.47490E-02 

4.41689E-02 

4.41716E-02 


l.OOOOOE+00 

2.61333E-02 

2.53455E-02 

2,53475E-02 



TABLE is; RAYLEGH PROBLEM Af=0.4 


y 

1st ORDER 

O 

1 

z 

3rd ORDER 

4th ORDER 

EXACT 

t»0.4 

.0 

l.OOOOOE+OO 

l.OOOOOE+00 

1 .OOOOOE+OO 

1.00000E*00 

1.00000E*00 

l.OOOOOE-Ol 

8.06701E-01 

8.U81BE-01 

8.10241E-01 

8,11062E-01 

8.11070E-01 

2.0C000E-01 

6.24995E-01 

6.33786E-01 

6.31A53E-01 

6.32468E-01 

6.32585E-01 

3.00000E-01 

4.64329E-01 

4.7448BE-01 

4.72i*88E-01 

4.73384E-01 

4.73290E-01 

4,OCOOOE-Ol 

3.30S85E-01 

3.39792E-01 

3.3BH33E-01 

3.38872E-01 

3.38980E-01 

5.00000E-01 

2.25617E-01 

2.32258E-01 

2.32390E-01 

2.32093E-01 

2.31998E-01 

6.00000E-01 

1.47779E-01 

1.51275E-01 

• 1.52066E-01 

1.51447E-01 

1.51494E-01 

7.00000E-01 

9.30980E-02 

9.37P.03E-02 

9.4b986E-D2 

9.43083E-02 

9.42643E-02 

8.00000E-01 

5.65815E-02 

5.53107E-02 

5.60055E-02 

5.58098E-02 

5.58295E-02 

9.00000E-01 

3.33039E-02 

3.10448E-02, 

3.14127E-02 

3.14501E-02 

3.14439E-02 

l.OOOOOE^OO 

1.90706E-02 

1 #660 1 6E-02 

1.67046E-02 

l#68l68E-02 

1.68274E-02 

II 






.0 

l.OOOOOE+00 

1 .OOOOOE + OO 

l.OOOOOE+00 


l.OOOOOE+00 

l.OOOOOE-01 

8.43748E-01 

8.A91B1E-01 

8.48E20E-01 


8.40767E-O1 

2.0000.0E-01 

6.93796E-01 

: 7.03-^46E-01 

7.02R54E-01 ' 


. 7.02917E-01 i 

3.00000E-01 

5.55<S27E-01 

5 »b81 4SF-01 

5.86947E-01 


5.67269E-01 

4.00000E-01 

4.33302E-01 

4.46431E-01 

4 • 45404E—0 1 


4.45601E-01 

5.00000E-01 

3.29126E-01 

3.40981E-01 

3.40286E-01 


3.40356E-01 , 

6.00000E-01 

2.43ft57E-01 

2.52890E-O1 

2.52365E-C1 


2,52559E-01 

7.00000E-01 

1.75983E-01 

1.R1954E-01 

1 ,R2022E-01 

UNSTABLE 

1.81926E-01 

8.00000E-01 

1.24163E-01 

1.26910E-01 

1 .27267E-01 


1.27124E-01 

9.00000E-01 

8.57025E-02 

; 6.57619E-02 

8.62730E-02 


8,61195E-02 

l*OOO00E*OO 

5.79671E-02 

i 5.61360E-02 

; 5.66541E-02 

. 

5.65305E-02 







TABLE 16! RAYLEIGH PROBLEM to=I.O, At =0.8 


I 


t 

y 

C -N 

3rd ORDER 

EXACT 

0.8 

,0 

l.OOOOOE+00 

l.OOOOOE+OO 

, l.OOOOOE+OO ■ 


l.OOOOOE-Ol 

6.36116E-01 

8.1612AE-01 

8.33029E-01 


2.00000E-01 

6.7R382E-01 

6.54357E-01 

6.73290E-01 


3.000COE-01 

5.32S10E-01 

5.183‘«9E-01 

5.27089E-01 


4.00000E-01 

A.03?6<5E-01 ■ 

4.01081E-01 

3.9907bE-01 


5.00000E-01 

2.93P35E-C1 

2.9H699E-01 

2.91841E-01 


6.00000E-01 

2.0587bE-0i 

2.12406E-01 

2.05903E-01 


7.000CCE-01 

■■ 1 .3RaH5E-0l 

1.439.3iE-01 

1.40016E-01 


8.00000E-01 . 

■ 8.951A6E-02 

9.30037E-02 

9.16903E-02 


9.00000E-01 

5.57005E-02 

5.7‘+]96E-02 

5.77798E-02 


l.OOOOOE+00, 

3.34646E-02 

3.39737E-02 

3.50152E-02 

2.4 

.0 

l.OOOOOE+00 

l.OOOOOE+00 

l.OOOOOE+00 


l.OOOOOE-01 

8.78H48E-01 

8.67418E-01 

8.78089E-01 


2.00C00E-01 

7.603B4E-01 

7.50ft77E-0l 

7.59006E-01 


3.00000E-01 

6.47J45E-01 

6.46011E-01 

6.48388E-01 


A.OOOOCiE-01 

5.41369E-01 

5.44P35E-01 

5.39498E-01 


5.00000E-01 

4.44847E-01 

4.45952E-01 ‘ 

4.43103E-01 


6.00000E-01 

. 3.5Rfll7E-0i . 

3.57136E-01 . 

3.57386E-01 


7.00000E-01 

2.8392^vL-01 

2.81284E-01 

2.82934E-01 


B.OOOOOE-01 

2.202bZE-01 

2.18A23E-01 

. 2.19768E-01 


9.00000E-01 

1.67402E-01 

1 .67061E-01 

1.67421E-01 


1 .oocooE>ao 

1.24599E-01 

1.25529E-G1 

1.25047E-01 




TABLE 17; 


RAYLEIGH 


PROBLEM At'O.I, t«24.9 


y 

C-N 

3rd ORDER 

EXACT 

0 • 

,10000E+01 

l.OOOOdE+OO 

.lOOOOE+Ol 

.lOOOOE+OO 

,95a76E+00 

9.54769E-01 

,95480E+00 

",200i)0£T00 

.90967E+00 

9.09682E-01 

,9097tlE + 00 

.30000E+00 

,.86i»86E + 00 

8.6488AE-01 

,86a97E+00 

.OOOOOE+OO 

,8P0a9E+00 

8.20516E-01 

,820baE+00 

.50000E+00 

.77668E+00 

7.76714E-01 

,7768t)E + 00 

".bOOOOt + OO" 

'.73357E + 00 

, 7.33609E-01 | 

.75379E+00 

,70000t+00 

.69126F+00 

6.91326E-01 1 

.69153E+00 

,80000E+.06‘ 

,6a993E+00 

6.49983E-01 

.65022E+00 

,90000E+Oa 

.60963E+00 

6.09686E-01 

,60995E+00 

.lOOOOE+Ol 

♦ 570il7F + 00 

5.70536E-01 

,8708aE+00 









no real advantage of the higher-order techniques. Further study may resolve 
some of the difficulties found here. 
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VI. RESUME- 


Polynomial spline interpolation has been used to develop a variety of 
higher-order collocation methods. Only those polynomials resulting in tri- 
diagonal, or at worst 3x3 block-tri di agonal , matrix systems have been evalu- 
ated. The governing systems are obtained directly in terms of the functional 
values and certain derivatives of the functional values at the specified nodal 
points. The development is generally given for a nonuniform mesh, for which 
a high degree of accuracy is maintained. 

Recently for a uniform mesh a so-called Fade, compact, Mehrstellung or 
Hermitian finite-difference procedure, which is 3x3 block-tridiagonal, has 
been proposed. ft is shown that this formulation is a hybrid method result- 
ing from two different polynomial splines. However, the Pade approximation 
is derived from a five-point discretization formulae and might be difficult 
to extend to nonuniform mesh systems. The hybrid spline results apply to 
variable meshes. Also, the compact system of equations does not include 
certain simple relations between the first and second derivative approxima- 
tions that are obtained from the polynomial spline interpolation formula. 

These relations are useful for reducing the size of the matrix system and 
thereby the computer time; in certain instances, boundary conditions can 
more easily be satisfied with these equations. 

Hermitian polynomial interpolation has been considered by Peters and it 
is shown herein that the compact differencing system Is derivable with this 
procedure. On the other hand, Peters is able to derive a tridiagonal matrix 
system involving only the functional nodal values Uy This would appear to 
be a significant improvement over the 3x3 compact system. However, the final 
system appears to be inconsistent, with an attendant loss of accuracy. 

Finally, from three-point Taylor series expansions and Hermitian discreti- 
zation of the functionals and their derivatives at the nodal points, an alter- 
nate derivation of the compact differencing scheme is presented. As only three 
nodal points are considered here, this procedure is less cumbersome than the 
Pade formulation and has been considered for nonuniform meshes and to develop 
systems with even higher-order truncation errors. Significantly, all of the 
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polynomial spline b lock- tri d 1 agonal systems can be recovered with this formula 
tlon. Moreover, a sixth-order (hybrid polynomial spline) 3x3 block- tri d i agona 1 
scheme has been devised. There does not appear to be any particular advantage 
of the polynomial spline formulation over the Hermitian discretization deriva- 
tion. 


The truncation errors of all the procedures are presented in tabular form, 

and results are shown for a variety of viscous flow problems- Of the fourth- 

f 7 81 

order methods, spline 4''* ' has the smallest truncation error. From the solu- 
tions to the model problems, the increase in accuracy with decrease in trunca- 
tion error is apparent. The sixth-order Hermi te formulation leads to extraordi 
nary accuracy even with very coarse grids. An important conclusion of the pre- 
sent study is that for equal accuracy the spline k procedure requires one- 
fourth as many points as a finite-difference calculation. This means less com- 
puter time and storage. 

Finally, polynomial interpolation has been used to develop higher-order 
implicit temporal integration schemes, which have previously been developed 
by Hermi te collocation. The present formulation leads to additional features 
of these methods, which are not obtained with the Hermi te procedures. For 
time-asymptotic (steady state) calculations there does not appear to be any 
advantage of these higher-order methods. For transient analyses, there are 
“stability" limitations not found with the lower-order techniques, and the 
integration is sensitive to small errors in the initial conditions. Therefore, 
for the examples considered herein, the advantages associated with higher- 
order time integration are minimal. 
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